Time-dependent density functional theory for quantum transport 
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Based on our earlier works [Phys. Rev. B 75, 195127 (2007) & J. Chem. Phys. 128, 234703 (2008)], 
we propose a formally exact and numerically convenient approach to simulate time-dependent quan- 
tum transport from first -principles. The proposed approach combines time-dependent density func- 
tional theory with quantum dissipation theory, and results in a useful tool for studying transient 
dynamics of electronic systems. Within the proposed exact theoretical framework, we construct a 
number of practical schemes for simulating realistic systems such as nanoscopic electronic devices. 
Computational cost of each scheme is analyzed, with the expected level of accuracy discussed. As a 
demonstration, a simulation based on the adiabatic wide-band limit approximation scheme is car- 
ried out to characterize the transient current response of a carbon nanotube based electronic device 
under time-dependent external voltages. 



I. INTRODUCTION 

The combined density functional theory (DFT) and 
nonequiHbrium Green's function (NEGF) approach has 
been widely employed to simulate steady state quantum 
electron transport through molecular junctions and other 
nanoscopic structures. In practice, DFT introduces an ef- 
fective single-electron reference system, on which NEGF 
analysis for steady current and electron occupation can 
be worked out. However, in principle it remains obscure 
how the conventional DFT for ground/equilibrium state 
of isolated systems can provide a rigorous framework for 
quantum transport, as electron transport is intrinsically 
a dynamic process. 

Time-dependent density functional theory (TDDFT)i 
has been developed to study quantum transport 
phenomenai^^— As a formally rigorous and numerically 
tractable approach, TDDFT promises real-time simu- 
lations on ultrafast electron transport through realistic 
electronic devices or structures. In early attempts, finite 
source-device-drain systems were treated by the con- 
ventional TDDFT for isolated systemsjii"— Since the 
source and drain were treated as finite, the results were 
not directly transferable to realistic devices coupled to 
bulk electrodes. By employing the NEGF approach and 
a partition-free scheme, Stefanucci and Almbladh have 
derived the exact equations of motion for the two-time 
Green's functions within TDDFT framework.^ They and 
their coworkers have further proposed a practical scheme 
in which the electronic wavefunction propagates in time 
domain subject to open boundaries,'^ The resulting nu- 
merical method has been tested on model systems. How- 
ever, its applicability to realistic electronic devices re- 
mains unexploitcd. 

Based on a reduced single-electron density matrix 
(RSDM) based formulation ) ^^d^ we have proposed a rig- 
orous TDDFT approach for open electronic systemsi^ 
We have also developed an accurate numerical scheme fS, 
based on a closed equation of motion (EOM) for the 
Kohn-Sham (KS) RSDM. Our RSDM based TDDFT- 



EOM has been combined with the Keldysh's non- 
equilibrium Green's Function (NEGF) formalism, the re- 
sulting TDDFT-NEGF-EOM approach has been applied 
successfully to exploit transient electronic dynamics in 
realistic molecular electronic devices.— One of recent ap- 
plications of the TDDFT-NEGF-EOM method was a 
simulation on the ultrafast transient current through a 
carbon nanotube based electronic device. It was found 
that the dynamic electronic response of the device can be 
mapped onto an equivalent classical electric circuit ^li^ 
which would be useful for future design of functional de- 
vices. 

Our previous simulations on realistic nanoelectronic 
devices^i^, have demonstrated the numerical feasibility of 
our TDDFT-NEGF-EOM approach for open systems, 
while its accuracy is largely determined by the quality 
of approximated exchange-correlation (XC) functional 
which accounts for the many-particle effects, and that 
of dissipation functional which characterizes the dissipa- 
tive interactions between the electronic device and elec- 
trodes. In our previous simulations we employed the adi- 
abatic local density approximation^^ (ALDA) for the XC 
functional and the adiabatic wide-band limit- (AWBL) 
approximation for the dissipation functional. Despite the 
success, the ALDA approximation is expected to become 
inadequate for addressing transient dynamics of open 
electronic systems in circumstances where electron cor- 
relations dominate. Moreover, the frequency dispersion 
part of the XC potential is completely missing from an 
adiabatic XC functional such as ALDA," which may lead 
to loss of crucial transient features in the electronic dy- 
namics. The AWBL approximation for dissipation func- 
tional may lead to nontrivial errors for electrodes of finite 
band widths and strongly inhomogeneous energy bands, 
or when non-Markovian memory effects play a significant 
role (such as in multi-channel conductors). 

Several groups have shown that within TDDFT frame- 
work the steady state current can be formally expressed 
by Landauer-Biittiker formula, if the steady state can be 
reachedi^i^ii^ This is important as the Landauer-Biittiker 



formula offers a convenient way to evaluate steady state 
current. However, in principle, it is the XC potential of 
TDDFT which explicitly includes frequency-dependent 
component that should be used, instead of that of ground 
state DFT. In practice, the use of an adiabatic XC func- 
tional would lose all the memory of transient dynamics, 
and results in the same steady state current predicted 
by ground state DFT. Therefore, it is desirable to go be- 
yond both the ALDA and AWBL approximations, and 
to develop more sophisticated XC and dissipation func- 
tionals for better characterization of transient electronic 
dynamics in quantum transport systems and structures. 

Time-dependent quantum transport has been ad- 
dressed from the perspective of open dissipative sys- 
tems. Rossi, Di Carlo and Lugli have developed a single- 
electron density matrix based Bloch equation to simu- 
late quantum-transport through mesoscopic devices J^ 
Burke, Car and Gebauer have proposed a single-electron 
density matrix based master equation for a current- 
carrying electronic system with dissipation to a phonon 
bathi? Weiss et al. have proposed an iterative path- 
integral approach)^ which is numerically exact but com- 
putationally expensive. Myohanen et al^ have ex- 
tended the Kadanoff-Baym approach to open interact- 
ing systems r^^ where the roles of initial correlations 
and memory effects were highlighted. The simulation 
of transient current has also been attempted by many- 
body quantum master equation approaches i^""— Yan 
and coworkers have started from a master equation 
based second-order quantum dissipation theory (QDT), 
and derived an EOM for the RSDM within TDDFT 
framework)^ where the bulk electrodes are treated as elec- 
tron reservoirs. This was followed by the recent estab- 
lishment of a formally exact QDT for the dynamics of an 
arbitrary non-Markovian dissipative systems interacting 
with bath surroundings.— The exact QDT is formulated 
in terms of hierarchical equations of motion (HEOM), 
and applied to solve time-dependent quantum trans- 
port problems .i^"-^ The HEOM-QDT is intrinsically a 
nonperturbative method, and is constructed to resolve 
the combined effects of the many-particle interaction, 
dissipative coupling strength, and memory time. The 
HEOM-QDT is by far the most tractable exact approach 
to time-dependent transient current through interacting 
electronic systems under arbitrary time-dependent volt- 
age. 

In this work, we aim at unifying our TDDFT-NEGF- 
EOM approach for open systems with the HEOM-QDT, 
to estabhsh a TDDFT-NEGF-HEOM method which 
is formally exact and numerically efficient. With the 
TDDFT-NEGF-HEOM approach, the dissipation func- 
tional readily goes beyond the AWBL approximation, 
resulting in systematic improvement on the simulated 
transient electronic dynamics. Moreover, the TDDFT- 
NEGF-HEOM approach is also applicable to devices op- 
erating at finite temperatures. 

This paper is organized as follows. In Sec.|ll] we 
briefly review the TDDFT-NEGF-EOM formalism for 
open electronic systems, and introduce the TDDFT- 
HEOM formalism. The two formalisms are then com- 
bined and resulted in a unified TDDFT-NEGF-HEOM 



formalism. In Sec. lIIIl we exploit formally exact compu- 
tational schemes to calculate the dissipation functional 
with the TDDFT-NEGF-HEOM formalism. In Sec.Hvl 
we propose three approximate schemes to reduce further 
the computational cost. In Sec.|V]a numerical example is 
presented, with analysis on the computational efficiency. 
We summarize our results in Sec.lVIl 



II. DENSITY MATRIX BASED 

TIME-DEPENDENT DENSITY FUNCTIONAL 

THEORY FOR OPEN SYSTEMS 



RSDM based TDDFT-EOM approach for 
isolated systems 



We start with ah initio many-particle Hamiltonian of 
an isolated system containing N electrons: 



H{t) 
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(1) 



where i and j run over all the N electrons, and v{ri, t) is 
the external potential for the i th electron due to nuclei 
and other electrostatic sources. Atomic units are adopted 
here and throughout the rest of manuscript. 

The Runge-Gross theorem^ establishes the one-to-one 
correspondence between the external potential v(r,t) and 
the time-dependent electron density function p{r, t) of an 
isolated system. A time-dependent Kohn-Sham (TDKS) 
scheme has been proposed to solve p{r,t), and hence all 
physical properties of the isolated system. The TDKS 
scheme follows the time evolution of an isolated reference 
system consisting of A'^ noninteracting electrons. The 
Hamiltonian of the reference system is as follows: 



N N 

H{t) = Y. ^(^*' ^) = E [ - 9^? + ^-ff (^*' *)] 



(2) 



where VcsifiTt) is the effective single-electron potential. 
A closed EOM has been derived for the KS RSDM a of 
the isolated system:--HL5 



ia{t)^[h{t),CT{t)]. 



(3) 



Here, h{t) is the KS Fock matrix. Its matrix element is 
evaluated via /i^t, (i) = j Xti-{'r)[~\'^'^+Vcs{r,t)]xv{r)dr 
in an atomic basis set {X/jI*")}- The square bracket on 
the right-hand side (RHS) of Eq. ([3]) denotes a commu- 
tator. The matrix element of cr is defined as (T^^(t) = 
(aj,(i)a^(t)), where a^(i) and aj^(i) are Heisenberg an- 
nihilation and creation operators for an electron occu- 
pying atomic orbitals n and v at time t, respectively. 
Fourier-transformed into frequency domain while con- 
sidering linear response only, Eq. ([3]) leads to the conven- 
tional Casida's equation^^ 



B. RSDM based TDDFT-NEGF-EOM approach 
for open systems 



In 2004, Fournais et al. have proven the real analyticity 
of electron density functions of any atomic or molecular 
eigenstates.^'' The real analyticity property leads to the 
holographic theorem of time-independent systems: The 
ground-state electron density on any finite subsystem de- 
termines completely the electronic properties of the entire 
system.— 

As for time-dependent systems, we have proven Holo- 
graphic Time-Dependent Electron Density Theorem'^ 
Let v{r,t) be the time-dependent external potential field 
on a finite physical system, p(r, to) be its electron density 
function at a given time to, and PD^f, t) the electron den- 
sity function within any finite subspace D. If p{r,to) is 
real analytic in r-space and v{r,t) real analytic in both 
r and i, there is a one-to-one correspondence between 
p]:,{r,t) and v{r,t); consequently, pjj{r,t) determines 
uniquely the full ab initio Hamiltonian, and thus, all elec- 
tronic properties of the entire time-dependent physical 
system. In principle, all one needs to know is the elec- 
tron density in a local subsystem. 

Based on the above theorem, we have developed a for- 
mally exact TDDFT-EOM formalism for open electronic 
systems. We consider coherent time-dependent quan- 
tum transport through a realistic electronic device. The 
entire system consists of bulk electrodes and the device 
region, where electron-electron scattering events mainly 
take place. For instance, in a two-terminal setup, the 
device region (D), left electrode (L), and right electrode 
(R) form the entire system. The region D is thus the 
open electronic system of primary interest. 

Expanded in an atomic orbital basis set, the RSDM 
(T of entire system can be partitioned into a number of 
blocks, among which cTa and ctq are the diagonal blocks 
corresponding to the electrode a (a = L or R) and the 
device region D, respectively; and (t^d = crl^^ represent- 
ing the off-diagonal block between the electrode a and 
the region D. The KS Fock matrix h can be partitioned 
in a similar fashion. 

To obtain the reduced electronic dynamics of the de- 
vice, we focus on the EOM of (Tq: 

icTD = [/iD,crD] -i^Qa(i). (4) 

a 

Here, Qa is the dissipation term due to the electrode a: 



Qa{t) — i {huaf^aD — (^Dahao) ■ 



(5) 



At first glance Eq. Q seems not closed. However, based 
on Holographic Time-Dependent Electron Density Theo- 
rem, all physical quantities are explicit or implicit func- 
tionals of the electron density of the system D, pD{r,t), 
and so is Qa{t). Therefore, Eq. ((4]) can be cast into a 
formally closed form: 

icr-D = [/iD [t;pT^{r,t)] , ctd] - i ^ Qq [t; pj^{r,t)] , (6) 



by noting that puifji) is the diagonal content of (TT:){t) 
in r-space. Equation ^ is the heart of RSDM based 
TDDFT-EOM approach. The transient electric current 
through the interface Sa (the cross section separating 
region D from the electrode a) can be evaluated via 



Mt) 



dr-^p{r,t) 



-tr[Q„(i)]. (7) 



The challenge now is to construct a formally exact and 
numerically accessible expression for Qa[i;pD(^,i)]-— 
Based on the Keldysh NEGF formalism, we have^ 



J — oo 

+ G<(i,r)SS(r,i)]+H.c. 



(8) 



Approximate schemes based on NEGF have been pro- 
posed in Ref. @. 

The exact HEOM-QDT provides another viable op- 
tion for the evaluation of Qa[t; pB{r,t)]. The result- 
ing TDDFT-HEOM approach is expected to improve 
systematically the accuracy of first-principles simulation 
on time-dependent quantum transport through realistic 
electronic devices. 



C. Hierarchical TDDFT for open electronic 
systems 

Differing from the conventional QDTs which assume 
weak device-electrode couphngs, the HEOM-QDT for- 
malism is considered as nonperturbative.^ We consider 
a reference noninteracting electronic system within the 
TDKS framework: 



Ht — Hd -\- 2_^{Ha + Hau)- 



(9) 



Here, Hd is the KS Hamiltonian matrix for the device 
D, Ha is the KS Hamiltonian matrix for the electrode a, 
and HaY) is the KS Hamiltonian matrix for the interac- 
tion between D and a. They are expressed as follows in 
second-quant iz at ion , 



-ffD = ^ h^^a^^a^, 



^iveD 



Ha — 2_^ ^ak d^f.dak. 



kea 



HaD — 2_^ 2^ ^akfidaf^lf^ + H.C. 
/^eD kea 



(10) 



Here, aj^ and a^, are creation and annihilation operators 
for the reference electrons in the device D, respectively. 
dl^f, and dak are the creation and annihilation operators 
for the reference electrons in the electrode a, respectively. 

V = {l^\kr,t)W) = J X^^{r)[-^^^ + v,s{r,t)]xAr)dr, 

eak = {ka\h{r,t)\ka), and takfi = {ka\h{r,t)\p). The 
EOM for the density matrix of the entire system (system 



plus electrodes), px, is 



Pt = -i[HT,pT] 



-«5m [-^"M^M + aU^M'PT], (11) 

with /„^ = J2k Kkf^dak and /t^ = J2k *«fc/^4fc- 

The TDDFT~EOM for the reduced system, Eq. ©, 
can be recovered. With the equality that cr^y{t) = 
trT[aJa^/5T(i)], Eg. pT|) leads to (c/. AppendixEl) 

id^ = [/iD, td] - Y. [V" W - V'^ W] , (12) 



where (pa,p,.(t) = trT[a+/Q^pT(^)] and v'l,pi/(*) = 

trT[/La^PT(i)]. 

By comparing Eqs. ([B]) and (|12p . it is apparent the 
dissipation functional Qa\pD{r,t)\ is directly associated 
with the first-tier auxiliary RSDM, (fait), as follows, 

= -iJde[<pc.{e,t)~<fiie,t)], (13) 

where (pa{t) ^ J de (^^(e, t) and ipl{t) = J de ip^e, t). 

A key feature of the exact HEOM-QDT is that for 
noninteracting systems, such as the TDKS reference sys- 
tem in the present case, the hierarchy terminates exactly 



J 



at the second tier without any approximation J^ The cor- 
responding HEOM for the KS RSDM and its auxiliary 
counterparts have been derived in Ref. [Ifl as follows. 



i(Paie,t) 



itpa,a'ie,e',t) 



+ [/„(e)-(TD]A„(e) 

+ X! / *Va,a'(e,e',0: 



(14) 



:-[e + A„(i)-e'-A„,(t)] ¥?„,«' 



(15) 

Here, /^(e) = l/[e'3(«-Mc) _)_ i] jg ^^e Fermi distribu- 
tion function for the electrode a, with /? being the in- 
verse temperature and pa the equilibrium Fermi energy; 
Aa{t) = —Va{t) is the energy shift for all single-electron 
levels in electrode a due to the time-dependent voltage 
applied on a; Aa,^^(e) = J2kea ^(<^ ~ e.k)t*„ktj^taku, is the 
device-electrode coupling matrix. 

Equations dm), dH]) and dH]) complete the TDDFT- 
HEOM formalism which is formally exact and closed. 
The basic variables involved are the reduced single- 
electron quantities {cr]3(i), y)Q,(e,i), tp^ ^./(e, e',i)}. 

D. TDDFT-NEGF-HEOM for quantum transport 

We now make connection between the TDDFT-HEOM 
formahsm to the RSDM based TDDFT-NEGF-EOM 
formalism. As shown in Appendix|Xl the auxiliary RS- 
DMs {ifia{e,t),(^a,a'{^,^' ^i)} cau be expressed in terms 
of NEGF quantities as follows^^ 



(^„(e,i)=W dr G<(t,r)S>(T,t;e)-G>(i,T)S<(r,i;e) , (16) 

J —OO 

(fia.,c.'{e,e',t)^i f dh f dt2\-Sc.'{t,h;e')GD{ti,t2)i:a{t2,t;e)Y (17) 

J —OO J —OO 

= 1 I dt, I di2{te(i,ii;e')G&(ii,t2) + S^.(i,ti;e')G<(ii,t2)ls>(i2,i;e) 

J —OO J — OO 

- S^,(i,ti;e')G^(ii,i2) + S>,(i,ii;e')GD(ii,i2)]s<(i2,t;e)}. (18) 
I 



Here, S|?;(i,T;e) are frequency-dispersed self-energies This is formally analogous to normal self-energies with 

{x = r,a,< and >), which normalize to normal self- S^ being the ground/cquilibrium-state quantities; see 

energies as / ^%{t, r; e) de = S^(i, r). It is straightfor- Eq. (|Bip . The following EOM thus applies: 
ward to derive from AppendixlB] that 

|sS(i,r;e) = -.[A„(i)+e]SS(i,T;e), 

^SS(t,r;e)= * [A„(r) + e] ^^(i, r; e). (20) 

Equations dUl and dSHUl) are the central equations of 
the TDDFT-NEGF-HEOM formalism. 



S;^(i,T;e) = exp 



i / A„(C)dC 



exp 



z / A„(C)dC 



e-''^(*-")i:S(e) 
tUt~T;e). (19) 



The stationary solution of TDDFT-NEGF-HEOM 
is related to the Landauer-Biittiker formula; sec Ap- 
pendix|B] for the details and remarks on this issue. 



III. HIERARCHICAL EQUATIONS OF MOTION 

SCHEME FOR CALCULATION OF DISSIPATION 

FUNCTIONAL 

A. Frequency dispersion scheme for 
TDDFT-HEOM 



In order to numerically solve the TDDFT-HEOM for 
Qa (i) , the integration over continuous function needs to 
be transformed into summation of discrete terms, i.e., 
f g{e) de — )■ X]a:=i Wkg{f-k)- Here, e^ and Wk are the /cth 
real frequency grid and its associated weight (0 < Wfc < 1 
-1 Wfe — 1), respectively. We have thus 



and Y.k= 



Nk 



Qa{t) = -i'Y^Wk ^ak{t) - ^ik{t) 



fc=l 



(21) 



Equations (|14p and P^ are recast into 



iiPak = [hnit) - efe - Aa(f)] ipak{t) + [fa{<^k) ^ Td] 
Wfc 
y. ^.a{tk)+^^Wk'iPak.a'k'{t), (22) 

a' k' = l 

i<fak,a'k' = - [Cfc + ^a{t) - ^k' " Aq/ (i)] iPak,a'k'{t) 

+ Aa'{ek')'Pak{t) - ¥3^,j.,(t)Aa(efc), (23) 

Hereafter, we adopt abbreviations ^Pakit) = fai^k^t) 
and ipak,a'k'it) = ipa,a'iek,ik',t). In practice, applica- 
tion of quadrature rules, such as Gauss-Legendre quadra- 
ture, reduces significantly the number of frequency grids 
compared to an equidistant sampling. This frequency- 
dispersed TDDFT-HEOM scheme applies to both zero 
and finite temperatures. Next we propose several effi- 
cient schemes particularly designed for finite temperature 
cases. 



B. Finite temperature schemes for 
TDDFT-NEGF-HEOM 



The ground/equilibrium-state self-energies S^(t) are 
usually exponentially decaying functions of time. There- 
fore, the following exponential expansion is attainable for 
t<>{T-t) a,tt>T: 



Nk 



tl{r-t)c,J2^'^ke^P[-lak{t~T)] 
fe=i 

^J2tl,ir-t) (24) 



k=l 



for X =< and >. jak are c-numbers with Ile{'jak) > 0. 
Moreover, the nonequilibrium self-energies S^j,(t, t) = 



e-i//A„(c)dc|^r!;^(^_^) Equation ^ becomes 



Nk 



fe=l 



(25) 



Here, ipak{t) is given by the RHS of Eg. (|T6|) with 
Sj!j(t, f;e) replaced by Sj^j.(r, t). The second-tier auxil- 
iary RSDMs, i^ak,a'k' (t) , are formally obtained by substi- 
tuting S^ {t2,t; e) with S^^ {t2 , t) on the RHS of Eq. ^ . 
EOM for ifak and ipak,a'k'it) now read: 

iifak = [huit) - ilak - Aa(i)] ipak[t) 
+ z[(TD(t)A>,^+a-D(t)A<fe] 



Nk 



+ X^2^ Vak,a'k'{t), 
a' k' = l 

iVak,a'k' = - [ilak + ^a{t) - ija'k' - Aq' (i)] <Pak,a'k 
A<^,)cpa,k{t) 



(26) 



* (^«'fe' 



^vlk'it)(.A>,~A<,) 



(27) 



Here, ctd = 1 — (Td is the KS reduced single-hole den- 
sity matrix. Eqs. (fT2|) . (|26l) and (|27|) constitute a refor- 
mulation of TDDFT-NEGF-HEOM formalism which is 
suitable for the numerical solution. 

Apparently, the computational cost for solving the 
TDDFT-NEGF-HEOM is determined by the number of 
exponential functions Nk used to expand the lesser and 
greater self-energies. Next we propose three decomposi- 
tion schemes to expand the self-energies as in Eq. ([24| . 



1. Matsuhara expansion decomposition scheme 

The ground/equilibrium-state lesser and greater self- 
energies can be evaluated via 



K{T-t) = ,^z /dee-(*-^)/^ne)Aa(6) 



(28) 



for t > T. Here, x =< or >, <r< = + and >;>=—, and 
/-(e) = fa{e) and /+(e) = 1 - /^e). Equation m 
states the fluctuation-dissipation theorem for electrode 
correlation functions. If the linewidth matrix Aa{£) can 
be approximated via a multi-Lorentzian expansion as fol- 
lows, 



A„(e) 



Na 

E 



Vd 



(e - f}rf)2 + Wj 



Aati 



(29) 



Here, rjd, ^d, and Wd > are the coefficient, center, 
and width of dth fitting Lorentzian function; and Aq^ is 
the corresponding frequency-independent linewidth ma- 
trix, respectively. The RHS of Eq. ([^5]) is calculated via 
complex contour integral and residue theorem, for which 
analytical continuation of Aa{e) and /^(e) into complex 
plane is needed. Aq(z) can be defined straightforwardly 
by replacing e with z on the RHS of Eq. pS)) , and f^ (e) is 
continued analytically to /^(z) in the same fashion. At a 



finite temperature, the Matsubara expansion scheme for 
the Fermi function is as follows, 



m^) 



1 






1 



■? 



Y^ 



^tt 



'^ap 



^ap 



(30) 



For the Matsubara expansion z^p = ^a + i7r(2p — l)//3. 

With all the poles of /^(z) and Aa{z) in the upper-half 
complex plane accounted for, we arrive at the following 
exponential expansion of self-energies: 



SS(r-t) 



with 






-lMt~T) 



Np 



-7cp(t— r) 



(31) 



ad 



^^^/-(^. + ^M/d)Ao 






/3 



lap 



-IZn 



-l/ia 



(2p - l)7r 



(32) 



We then make connections between Eq. ([31]) and Eq. ([24| . 
The total number of exponential terms is the sum of 
number of Lorentzian functions and number of Matsub- 
ara terms considered, i.e., f^k = ^d + -^p- In principle, 
7Vj3 — >■ +00 is required to achieve an exact expansion for 
S^'>(r — i). In practice, a smallest possible Np but guar- 
anteeing an expected accuracy is desired. The value of 
Np increases rapidly as the temperature is lowered. 



2. Partial fractional decomposition scheme 

A major disadvantage of the Matsubara expansion is 
the poor convergence at low temperature. Alternatively 
we may adopt the partial fractional decomposition (FED) 
of Fermi function. "^^ The FED expansion is formally iden- 
tical to Eq. (I5n)) . but with Zap = Ha + 2y^//3. Here, Cp 
is the p th eigenvalue of the Np x Np matrix Z defined as 
follows. 



2m{2m — 1) 5„ 



-2Np{2Np-l)5^.,N,. (33) 



with the constraint that Im(y/e^) > 0. 

With the same multi-Lorentzian expansion of Eq. ([29| , 
the self-energies are decomposed into exponential func- 
tions exactly as Eq. ((3T|) . The FED scheme leads to a for- 
mally similar exponential expansion for the self-energies, 
and hence for the resulting HEOM, as compared with 
the widely used Matsubara expansion scheme. However, 
numerical tests have confirmed that with same Np, the 
FED scheme yields much more accurate approximation 
for fa{z) than the Matsubara expansion scheme. '^^'■^^ In 
other words, to achieve the same level of accuracy for 



the self-energies, a much smaller number of exponential 
terms is required with the FED scheme. Therefore, the 
FED scheme is superior to the conventional Matsubara 
expansion scheme in terms of computational efficiency. 



Hybrid spectral decomposition and frequency dispersion 
scheme 



As mentioned after Eq. (1321) , in principle Np — )■ -t-oo 
is required to achieve an exact expansion for the Fermi 
function, and hence for the self-energies. With a finite 
Np, the deviation between the RHS of Eq. (15(11) and the 
exact function /^ (z) results in the following contribution 
to the self-energies: 



■ix.dcv 



(r-i) 



^xl 



Im(z)=y 



e*^(*-"V^i^)A„(z)dz (34) 



for i > T, as shown in Eig.[T] Different from Eq. 
where the integration is along the real axis, here the in- 
tegration contour is a horizontal line: Im(z) = y > 0. 
The value of y is chosen so that the first upper-plane Np 
Matsubara poles reside in the interstitial region bounded 
by Im(z) = and Im(z) = y, while all the other poles 
are located outside this area, i.e., {2Np — l)7r//? < y < 
{2Np + 1)^//?. 




Re(z) 



FIG. 1: An illustrative plot showing the contour of integra- 
tion for fa{z)A.a{z) with a single-Lorentzian spectral density 
of fii = 5. The dots equidistantly distributing on Re(2) = 
correspond to the Matsubara poles (grey dots) in the up- 
per complex plane. Integration along the dark solid arrow 
lying right on the real axis gives the desired self-energy of 
Eq. pSfl . However, the resulting HEOM encounters numeri- 
cal difficulties. Due to the fact that the two dark solid arrows 
form a closed loop (completed with the dark dashed arrows at 
infinitely distance carrying zero values), the self-energy can 
be obtained alternatively by summing up the residues at the 
poles within this loop. 

As long as the modulus of Aq(z) remains sufficiently 
small along the line Im(z) = y, Eq. (p4|) presents a com- 
plementary contribution to the RHS of Eq. ([31]) . The 
integration in Eq. p4p is transformed into summation 



of discrete terms, following the frequency dispersion 
scheme introduced in Sec. lIII Al i.e., L , x g{z)dz — > 

J2q=i ^g 9i^q + *y)- We have thus 



N„ 



SS^<i-(r-i)^^i^S^e-^-(*--), (35) 

9=1 



with 



K^q = <^xiWqfl^{eq + iy)Aa{€q + iy) , 



(36) 



Therefore, overall S^(t — t) is given by combining the 
RHS of Eqs. ^ and ^. The total number of expo- 
nential functions is Nk = Nd + Np + Nq. In practice. 



to make use of Eq. ([35|) . one needs to assign an appro- 
priate number of Matsubara terms, Np. On one hand, a 
smaller Np leads to fewer unknown variables of HEOM. 
On the other hand, a too small Np (and hence a too 
small y) would cause severe convergence problem when 
solving the stationary solutions of HEOM. Therefore, the 
value of Np should be chosen carefully for the balance be- 
tween computational cost and numerical difficulty. Once 
the Np and y are settled, efficient quadratures can then 
be adopted to optimize the frequency dispersion (the set 
{wq, eq} with q = l,...,Nq) for Eq. (jMl)- Finally, all the 
Nd poles of Aa{z) inside the region < Im(2:) < y are 
collected. 



It is important to note the value of Nd depends on the 
mathematical form of basis functions, based on which 
A„(e) is expanded. The widely used multi-Lorentzian 
expansion has been adopted in Sec. lIII B"T] see Eq. (P^. 
Each Lorcntzian function gives a pole in the upper com- 
plex plane. One of the advantages of using a Lorentzian 
function ghi^) is that the magnitude of its analytically 
continued counterpart, |5l(^)|, decays smoothly as Im(z) 
increases. Therefore, the deviation term, S^''^''^(r — t) 
of Eq. ([34| , becomes consistently less significant as more 
Matsubara poles are considered explicitly. Alternative 
basis functions can also be used. For instance, a Gaus- 
sian function is analytic, and does not have any pole on 
the entire complex plane. Therefore, for Aa{e) fitted by 
a number of Gaussian functions, Nd — 0, which reduces 
the number of unknown variables of HEOM. However, 
care must be taken when treating the "deviation" con- 
tribution, as some of the K^ may assume large values, 
due to the nontrivial value of complex Gaussian function 
at certain z — eq + iy. 



IV. APPROXIMATE SCHEMES FOR 
CALCULATION OF DISSIPATION FUNCTIONAL 



A. Schemes based on wide— band limit 
approximation 



Even with the sophisticated hybrid scheme for decom- 
position of self-energies, the exact HEOM approach for 
the reduced electronic dynamics can still be expensive. 
Efficient approximate schemes are thus required. For 
clarity, we will omit the spin index in the following deriva- 
tions of approximate schemes. 



1. TDDFT-NEGF-HEOM formalism under wide-band 
limit approximation 



The wide-band limit (WBL) approximation involves 
the following assumptions for the electrodes: (i) their 
bandwidths are assumed to be infinitely large; and (ii) 
their linewidths are assumed to be energy-independent, 
i.e., A.a{e) « A^. Denote A^ = ttAq. From the def- 
inition of self-energies and the expansion of the Fermi 
function [Eq. (PU)) with Np — )• -l-oo], one obtains for t > r 

S<'>(T,i)^±*e'/^''^'^^°(*^)y'd6/±(e)e"(*--)A„(6) 



±.i5{t — t)Aq, 
2 



±^e»/:rf*i-=.(*i)A„, 



P 



(37) 



k=l 



where Zak{t) = Zak + ^a(t)- The expansion of Fermi 
function thus leads to a sum of exponentials for the self- 
energies. We introduce auxiliary self-energies Safe to 
account for the exponentials, i.e., 

S<^> (r, t) = ±iS{t - t)A„ + Y, ^o^kir, t), (38) 



Safe(r,i) = -e'^.*'**i"-(*i)A„ 



(39) 



It is implied that 'Eak{t,t^) — IAq. Next, we insert 
relation Eq. (l38l) into Eq. (IT6)) . and arrive at 



¥>a {t) =l[l-2(Tu{t)]Aa+Y, Vo^k {t) , (40) 



Numerical tests on model systems^Si have confirmed 
that as the temperature lowers, the total number of ex- 
ponential functions needed to accurately resolve the self- 
energies in the hybrid scheme grows much slower than 
that in the Matsubara expansion scheme. 



with the auxiliary current matrices being as follows, 
<^afc {t) = -i I dT [G^ (i, r) - G^ (i, r)] Safe (t, t) 



dTGl,{t,T)'S^k{T,t). 



(41) 
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Under the WBL approximation, NEGF formalism gives 



G^j,it,T) = -^^{t - t) W^{t) W+{t) 



W^{t) = exp± < =Fz / dti hoiti) - iA 



(42) 
(43) 



Here, we take the time from which the voltage emerges as 
to = 0, i.e., Aait < 0) = 0, and A = J2a^o,- Therefore, 
the EOM of tpak{t) can be established readily as follows. 






huit) - iA - Zc,k{t) 0ak{t). (44) 



The coupled equations of motion, Eqs. (jj]), (l25t . (|40)) and 
(|33]) can be solved together for a complete description 
of the nonequilibrium electronic dynamics of the device. 
Equation (|M1) is self-closed, which suggests that under 
the WBL approximation the TDDFT-HEOM automat- 
ically terminates at first tier, instead of second-tier ter- 
mination for the exact formalism. Thus the additional 
second-tier auxiliary matricesiS are not needed with the 
WBL. Similar apporach has been adopted by Croy and 
Saalmann<^ 



integral within the square bracket on the RHS. The basic 
idea is to first disregard the time-dependence of ft.D(0 
and Aa{t) in Eq. (|50|) . evaluate the time integral, and 
then restore subsequently the time-dependence of Hd (t) 
and Aa{t). The final approximated expression for Pa{t) 
is as follows. 



P„(<)^— C/+(i) 

TT 



-!-oo 



&/a(e)e* 



1 



e-hB{0) + iA e- ft,D(i) + JA + Aa(t) 



CXD 



e-hr,{t)+iA + Aait) 



de}Aa. (51) 



The RHS of Eq. dCT]) gives exact Pa (i) for the ground 
state at i < and the steady state at t ^ +oo, provided 
that Aa{t — >■ -f(X)) is a constant. At < t < -l-oo, it 
provides an adiabatic connection between the ground and 
steady states. Equation ([5T|) can be solved conveniently, 
along with the EOM for U+{t) as follows, 

I U+ (t) = [h^ (t) - ^ A - A„ (i)] U+ (t) . (52) 



2. TDDFT-NEGF~EOM formalism under adiabatic 
wide-band limit approximation 

We consider another approximate scheme for Qa(t) 
based on WBL approximation for electrodes, as well as an 
adiabatic approximation for memory effect. The scheme 
aims at simplifying the TDDFT-NEGF formulation for 
Qa {t) ; see Eq. ^ . Due to the WBL approximation, the 
ground/cquilibrium-state self-energies become simply 



-T)Aa, 






(45) 


.J —oo 


-^Ue 


A„. 


(46) 



S<(T,i)- 

By inserting Eqs. p5|) and (|46| into Eq. ([8]), the dissipa- 
tion functional is formally simplified to be 

QAWBL(^) = {A„,<Td} +Po(t) + [Pa{t)]^ , (47) 

/+00 
G[,(i,T)S<(r,i)dT. (48) 

-oo 

Before the external voltage is applied, the entire compos- 
ite system is in its ground/equilibrium state. We have 



Pa(i) 



■u+it) 



+ 00 



/o(e) 



1 



lie 
la 



.. — 4€'?"rr — 



U-{r)dT 



e-huiO) +iA 
}a„. (49) 



e"*de 



t/±(t)=cxp±|Ti/ dr[/iD(r)-iA-A„(r)] } 



(50) 

To facilitate the calculation of Pa (t) given by Eq. (j49| , 
an adiabatic approximation is adopt to evaluate the time 



B. Scheme based on complete second— order 
quantum dissipation theory 

As presented in SecUHl the TDDFT-NEGF-HEOM 
for KS RSDM and associated auxiliary matrices termi- 
nate exactly at the second tier. This amounts to an ex- 
plicit treatment of device-electrode interaction at leading 
4th-orderiiS The number of unknown matrices at zeroth, 
first, and second tier is 1, NaNk, and N^N^, respec- 
tively. Here, N^ is the number of exponential functions 
used to resolve the memory contents of self-energies or 
electrode correlation functions. Obviously, the second- 
tier variables dominate the computational cost for solv- 
ing the TDDFT-NEGF-HEOM. Therefore, a straight- 
forward way to reduce the computational expense is to 
avoid explicit involvement of the second-tier variables, 
i.e., to truncate the TDDFT-NEGF-HEOM at first tier 
approximately. 

A simple truncation scheme is to set all second-tier 
auxiliary matrices to zero, i.e., ^ak,a'k'{t) = 0. This 
corresponds to the chronological ordering prescription of 
2nd-order QDT. 

An alternative 2nd-order approximation is as follows: 



G<{t,r) 

G>{t,T) 



l<Tj,{t)U+(t,T), 

-tau{t)U+{t,T), 



with 



U^i't^'r) =exp± 



T^ 



hMQdC 



(53) 



(54) 



The dissipation functional becomes approximately 

Qa (t) « [cTD (i) n> (t) + &u (t) n< (t) + H.C.] , (55) 



where 



Ul(t)^lf dTU+(t,T)i:l{T,t). (56) 

Following Sec. lIIIi we decompose self-energies as linear 
combinations of exponential functions; see Eq. (|24|) . This 
leads to 

ns(t) = 5^n^,(t), 
fc=i 



J —oo 

The EOM for n^fe(i) is thus 



(57) 



tUl,{t) = [huit) ~ Z7afc - A^{t)]Ul,it) - Al„ (58) 



with A^j, and 70,^ referring to Eq. (|24|) . 

The resulting complete 2nd-order QDT approach thus 
involves the coupled EOM for (TY){t) and {Il^^(t)}, 
with the total number of unknown matrices being N = 
2NaNk + l, much fewer than N^N'^ alone. This approach 
can be improved further, by formal inclusion of higher- 
order device-electrode interaction into the reduced sys- 
tem propagator U^{t,T). This can be achieved by at- 
taching self-energy to hj:,{t) on the RHS of Eq. ([54|) . 



V. NUMERICAL RESULTS 

As a test for our practical schemes, we have performed 
the numerical calculations employing the AWBL approx- 
imation as outhned in Sec. IV A. 2. LODESTAR^^ was 
used to carried out the calculations. The system of inter- 
est is a (5,5) carbon nanotube which is covalently bonded 
between two aluminum electrodes, as shown in Fig. [21 In 
the simulation box, we include a finite carbon nanotube 
and 32 aluminum atoms for each electrode explicitly. The 
minimum basis set ST0-3G is adopted in the calcula- 
tions. We simulate the time-dependent electric current 
through the left (right) electrode by taking the trace of 
the corresponding dissipative term Ql (Qr); see Eq. ([T]). 

Figure |3] shows the currents versus time for different 
numbers of carbon atoms, i.e. 20, 40, 60 and 80. The 
bias voltage Vb is switched on exponentially a.t t = 0, 
and as shown in Fig. [31 We observe that the currents 
reach their steady states in about 8 to 12 fs for different 
systems. Also, the value of the steady state current does 
not increase proportionally with the size of the carbon 
nanotube, as that of classical case. We note from Fig.[3l 
that the steady state currents for the systems with 20, 
40, 60 and 80 carbon atoms are about 21, 9, 14 and 4 
^A, respectively. This shows that the quantum finite size 
effect plays an important role at such small scales. We 
point out that all the above calculations are performed on 
single-CPU desktop personal computer, and the memory 
is 2 gigabytes. 




FIG. 2: The ball-and-stick representation of the system of in- 
terest which is a carbon nanotube (5,5) welded to aluminium 
electrodes. There are 60 carbon atoms for the carbon nan- 
otube in this case. 



VI. CONCLUDING REMARKS 

We have proposed a first-principles TDDFT-NEGF- 
HEOM method for transient quantum transport through 
realistic electronic devices and structures. The TDDFT- 
NEGF-HEOM formalism is in principles exact, and for- 
mally equivalent to the TDDFT-NEGF-EOM formalism 
proposed previously. In practice, it involves hierarchical 
equations for KS RSDM of reduced system and associ- 
ated auxiliary quantities. 

By resolving memory contents of self-energies or elec- 
trode correlation functions, we construct QDT-HEOM 
in the KS reduced single-electron space. The resulting 
TDDFT-HEOM exactly terminate at the second-tier. 
TDDFT-NEGF-HEOM formalism combines TDDFT- 
NEGF-EOM and TDDFT-HEOM, and its computa- 
tional cost is determined by the number of exponential 
functions used to expand the self-energies. Various de- 
composition schemes are presented. 

To reduce computational cost further, we have also de- 
vised approximate schemes for TDDFT-NEGF-HEOM. 
One is based on the WBL treatment for self-energies. 
Another one is based on complete 2nd-order QDT. These 
approximate schemes significantly reduce the computa- 
tional cost, and improve the efficiency for solving the 
reduced electronic dynamics. Additional adiabatic ap- 
proximation is introduced for the WBL scheme, and the 
resulting AWBL approximation turns out the most effi- 
cient scheme. Numerical simulations based on the AWBL 
scheme demonstrate that first-principles simulation can 
indeed be carried out for real devices, and the interesting 
results have been obtained. 
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Appendix A: Detailed derivation of HEOM for an 
open noninteracting system 

Here we give a direct derivation of HEOM for a nonin- 
teracting system, without resorting to the path-integral 
formalism for Fermion operators. Similarly with Ref. lid . 
we introduce the reservoir _ffB~interaction picture, where 
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Time (fs) 

FIG. 3: The grey solid line represents the bias voltage applied 
on the systems. The bias voltage is switched on exponentially, 
Vb = Vo(l - e"*/°) with Vo = 0.1 V and time constant a = 1 
fs. The transient currents for the systems with 20, 40, 60, 80 
carbon atoms are shown dark dotted line, dashed line, solid 
line and dash-dot-dot line respectively. 



Hb = X^a ^a ■ In this interaction picture, the total 
Hamiltonian is 



HT{t) = HY, + Y,HMt)- 



(Al) 



Here, -ffcoW = E^ [/Im(^)«a' + o-lMt)] where 
flfi (0 = J2k^aktj.dak{t)- The quantum Liouville equa- 



tion reads 



= -i/:DPT(i)-J^/:aDPT(t), (A2) 



where the superoperators are defined as 

CaD ■ = [HaDit), -J. 



(A3) 



First, we establish the EOM for the RSDM of primary in- 
terest, crDp,y(t) = tr^[aJ,a^pT(i)] = tr^[aJ,a^pT(0]- Fol- 



lowing Eq. (|A2p . we have 



a 
a 

Ai'eD 
^ [^D, TD(i)]^, - E [v=c.Mt) - 'pI^M ■ (A4) 



Here, the second equality results from trace cycling in- 
variance property; and the fourth follows the definition 
of auxiliary RSDM, ipa,f^uit) = ti:^[fl^{t) a^pT{t)], and 
its Hermitian conjugate <pj,^^,,(i) = tT^[alfaf,{t)pT{t)]. 

Considering the multiple-frequency-dispersed scheme 
in Ref. llCt we can introduce the frequency dependence 
of / and /' from the beginning. In the reservoir Hb- 
interaction picture. 



dUt)^-¥c.kitlHB] 

Simple integration leads to 

dL(t) = e^-^4[-+^=Ml-Jt^, 



(A5) 



k£a 






X e 



./4drA„(r) 



(A6) 



Let 



fU^. i) = ^ <5 (c. - £„,) e-(*-*") t^k, dl 



k£a 



X e 



z/4drA„(r) 



(A7) 



we ha,ve fl^{t) = J duifl^{uj,t), /„^(t) =J dujU^{u,t), 
and HaD{t) = E^ / duj[flf,iuj, t)af, + a^Jat.{^, t)] . 

Define 

^a,^u{i^,t) = tv^{fl^{uj,t) a^,pT{t)] , (A8) 

^l,t.A^,t)=ir^{alfo,^{uj,t)p^{t)]. (A9) 

They satisfy ipa^^v{t) == J duj ipa^f,^{uj,t) and ipi^^^{t) = 
/ duj (fl^ u,^{^, t). Comparing with the NEGF formalism. 
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/oo 



it is readily verified that where the second equality follows the derivation for 

Eq. (B7) in Appendix B of Ref. IH. With the Langreth's 
(Pa^^i,{u}, t) = —i \J 5 {uj — Eak) G'^ak (^j ^Qfcjy- (AlO) analytic continuation rules,— we have 

The following nonequilibrium Green's function is defined 
on the Keldysh contour, where the time variable goes 

from — oo to +00 then back to — cx): 

G^^k {t,t') = -ztr, {TcK M4/C ir')]pT (to)} 

= 53 / '^'^1<^M>' ('T' 'ri)t*aku 9ak (ti, r') , 

(All) 



+ G<,(i,ti)5°Jti,t)]. 
Inserting Eq. (|A12p into Eq. (jAlOp leads to 



(A12) 



/OO 
dti [Gl,,, {t,h)g<,{h,t) + G<^, (t, ti ) g^fc (ti , i)] 

/oo 
dti [G;,, (t, ii) S< ,,, (ii, i; 0.) + G<,, {t, h) Y,ly, {h,t; ^)] 

^iY. I dh[G<^,it,h)i:l^.^{h,t;u;)~G>^,it,h)-E<^^,^{tut;u;)], (A13) 



i/'eD" 



where the third equality employs the relation 

X'-^'^ (t,r) = ±1? [± (i - r)] [X>(i,r) - X<{t,T)] 

(A14) 
with JC being Gd or S^. The frequency dispersed self- 
energies are defined by 



kGa 

(A15) 
where x = r, a, <, and >. They satisfy 

= i 22 Kiku'^akv / duj 6{UJ - eak)fa {^ak) 
kea •' 

= ^ /"dL^Aa,.v(w)/a(c^)e^"(*-*i) 

xe'^^'^'''^-^'\ (A16) 



RSDM, ipa,f_iiy{uj,t), as follows, 

= tl'T {/L('^> t) a^ PT(t) + /L(w, t)a^ PT(t)} 
= i[uj + Aa{t)] ipa,^i.iuJ, t) 

+ trT{/L(^> t) a^ [ - iCr, PT{t)] } 

+ Yl ^^T{fL{^^ t) a^i [ - iCa'D PT{t)] } 

a' 
= i[uj + Aa{t)] ipa,f^i,{uJ, t) 

- « tr.j, { [fl^iuj, t) a^, Hu] PT{t) } 

- f ^ tr^ { [fl^{uj, t) Qf,, Ha'D{t)] PT (i) } 

a' 
^i[uj + Aa{t)] ipa,f^,y{uj, t) 

-i Y hD,j.,^'tT^{fl^{uj,t)af,, pT{t)} 
M'eD 

a' 

+ i^ Y trT{*afc,.**fcp''5(a;-eafc)aj^, a^pT(t)} 

kea fj,'eD 

= i[uj + Aa(t)] ipa,f^u{(^, t) 

/i'GD 

-zYldio'tT,{fl{io,t)U^{u;\t)pT{t)}. (A17) 



Next, we establish the EOM for the first-tier auxiliary 
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Comparing with the NEGF formahsm, it is readily ver- 
ified that 

= ^ ^ S{UJ - eak) S{uj' - ea'k')taki^tl,f.,^ 
k^a. k' Ga' 

= -i ^ ^ 5{iO- Cak) S{uj' - ea'k')takvt*a'k',i 
k(Ea k' (Ea' 



Inserting the second term on RHS of Eq. (IA19P into 
Eq. (JAISP while using Eq. (|A15|) . we arrive at a term on 
the RHS of Eq. (IATtI 



<Paa',iyfii^,^',t) = -l 



V" d,Ti dT2['Ea',f^f^At,Ti;Uj') 

.. ^T^JC Jc 



l-Lll.L-2£'D 



X G^i^2(Ti,r2)Sa,^2^(T2,t;a;)] 



(A21) 



xG<,,„,(i,i)- 



(A18) 



Generally, the NEGF on the same Keldysh contour sat- 
isfies 

Ga'k\ak {t, t') = -i tr^ {Tc [da'k' (t) c^lfc(T')] Pt(^o) } 

— Saa'Skk' gak{T,T') 



+ y~] / dri / dT2ga'k' {T,Ti)ta'k'fj.i 
.. .. ^T^Jc Jc 



tiifj,2eD 



It turns out that the EOM for the first-tier auxiliary 
RSDM, ipa,fj,,y{(^,t), is formally identical to that given 
by path-integral formulation in Ref. lid . Combining 
Eqs. (|JT8l) through ([MB . Eq. (Kl7\ finally reads: 

i^a,^iy{'-0,t) = - [w + AQ(i)]y>a^p^(w,t) 



X Gpips (ri,T2) t^s,^^ 5afc (7-2,1-') , 

(A19) 

where the same technique in Ref. [S^ adopted for deriva- 
tion of Eq. (lAlip is used. At r = r' = i the lesser 
component of Green's function for the uncoupled lead 
is g<^ (i,i) = g<^{t - t) ^ if a {tqk ). Ins erting the first 
term on RHS of Eq. ((M9)) into Eq. (|JT8)) leads to 

- « ^ ^ ^ I dio' 5{lO- Cak) S{uj' - Ca'k') 
a' k^a k' ^a' 

X *aku*^a.'k'ii^Cia'^kk' \}ja \^ak) \ 

= Aq^^,{uj)fa{uj). (A20) 



J2 duj' iPaa',^^,{uJ, cj\ t). (A22) 



Here the compact definition ()A2ip establishes the natural 
connection between the second-tier auxiliary RSDM and 
NEGF quantities. Again, using the Langreth's analytic 
continuation rules,— we can immediately write down a 
real time expression for Eq. (JA21I) as follows, which looks 
however more complicated. 



^ \ dt2 dTiSa',^pi(t,Ti;w')G^iP2(Ti,i2) 

/ dt2 / dTi'Ea',fj.fj.i{t,Ti;Lu')Gfj.^f^^{Ti,t2) 

J-oo Uc 



A'lA'2eD 



S> .(t2,i;a.) 



ic 

nt nt 



S< ,(t2,t;^) 



= i J2 { [ ^*2 / dii[S<^^^^(i,ti;w')G;:,^,(ii,t2) + S;',^^,(t,ii;^')G<,,,(ii,t2)] S>^,,(t2,i;^) 

f dt2 I dti[S^,,^^,(t,ii;^')G>^,(ii,i2) + S>,,^^,(i,ii;^')G;:,,,,(ti,i2)] S<^,.(i2,i;^)j. (A23) 

J —oo J —oo ) 



M1M26D 



The last equality recovers Eqs. (22) and (24) in Ref. |40|. 
At t ^- — cx) the device and leads are fully decoupled, i.e. 



Pt(— 00) = pB <Xi pd(~oo). Noting that 
AQ,^^(w)/a(a;) 

duj' tv^ [fLi^, t)fa'p.{uj\ t) pt(-oo)] , 



(A24) 
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we can have a most simple and compact expression of the Based on Eq. (jA25[) . it is easier to obtain the EOM for the 
second-tier auxihary RSDM as follows, second-tier auxiliary RSDM, faa' ,u^.{'^: ^' , t), as follows. 



= tT^{fU^,t)U,^{Lu',t) [ptW - Pt(-c^)] }. (A25) 



¥'aa',i/^(a;,a;',i) 

== i[uj + Aa{t)] ipaa',iyfi{l^,l^',t) - i [uj' + Aa' (t)] (faa'^iyfiiuJ , Uj' ,t) + tT^[fl^{uJ,t) fa' f^{uj' ,t) pT{t)] 
= i[uj + Aa{t) -Uj' - Aa'{t)](Paa',uf,{uJ,Uj',t) + tT^{fl^{uj ,t) fa' f,{uj' ,t) [-iCnPTit]] } 

ai 
= i[uj + Aa{t) -Uj' - Aa'{t)] ipaa',v^i{^,^\t) - i tr^ { [/^^ (w, i) /q'^ (w', i), iJo] PtC*) } 

- iY,tr^{[fi,{uj,t) fa'^,{uJ' ,t) ,Ha'i,{t)\ -pTit)] 

Oil 

= i[u> + Aa{t) ~Uj' ~ Aa'{t)] ipaa'.^fj^iuJ,Uj',t) 

-«X!X! Yl Yl ^^T[fL{^,'t)0'f,'Sa'aiSk'kiS{uj' ~ ea'k')taiki,,'t*a,k'f,P'T{t)] 

+ ^J2Y1 Y XI *^T [a-l'U^iu}': t) 6aaiSkkiS{uj - eak) taikiu t*ak^' PT{t)] 

Oil /.i'GD A,'i Gcti k^a 

= i[uj + Aa{t) -LU' - Aa'{t)] Ifaa' ,i^f^{^ , ^' ,t) -i ^ Aa' ,fj.fj.' {ui')ipa,fj.' „{lO ,t) +i ^ If]^, ^^^, {uj' ,t)Aa,fj.'i,{Lo) . 

ij,'eD ii'eD 

(A26) 



The last equality recovers Eq. (6.6) in Ref. [Tol 



Appendix B: Remarks on Landauer— Biittiker 
formula for steady current 



Under an external voltage applied to electrode a, a 
homogeneous energy shift A^ is resulted for all energy 
levels, and the self-energy is thus 



S^(i,T)=exp 



AaiOdC 



^Ut-r), (Bl) 



where x = r,a,<, and >; and the tilde symbol denotes 
quantities of ground/equilibrium-state in absence of any 
voltage, where translational invariance in time applies. 

For the composite system to approach a steady state 
as i — >■ +00, it is necessary to have A^ = Aa{t — >■ +00) 
approaching a constant for any a. Based on Riemann- 
Lebesguc lemma, Gd(^, t) — Sa(i,T) = as i — > +00 
and r remains finite. For t,T ^ +00, Gu{t, t) = GD{t — 
r) and SQ(i, r) = Sq,(< — r), i.e., translation invariance 
in time is retrieved at steady state. 

The NEGF formalism gives steady current through 



I 

electrode a, J^ = Ja{t ~> +00), as follows, 

= E Jde[fa{e)-fa'{e)]Taa'{e), (B2) 

Taa'ie) = 27rtrD[c?^(e)A„,(e^)GS(e)A„(e^)]. (B3) 

Here, e^ = e — A^, and Taa'i^) is the KS transmission 
coefficient between electrodes a and a' . Equations (|B2p 
and (IBSp correspond formally to the Landauer-Biittiker 
formula. Equation (jB2[) has also been reached in the 
framework of HEOM-QDT, see Ref. ^ It is worth em- 
phasizing that the exact XC potential of TDDFT is in 
principle frequency-dependent. 

Equations (jB2[) and (jB3p have been widely employed 
in the DFT-NEGF approach to simulate steady current 
through molecular devices and structures. The exact XC 
potential may depend explicitly on the time evolution his- 
tory of electron density in TDDFT, while it is determined 
exclusively by p{r,t — > -l-oo) in stationary DFT. Un- 
der the same time-dependent applied voltage, DFT and 
TDDFT frameworks result in explicitly different steady 
current if and only if a nonadiabatic (or frequency- 
dependent) XC potential is adopted in TDDFT frame- 
work. In other words, DFT-NEGF approach is regarded 
as an adiabatic version of TDDFT-NEGF. 
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In principle, the steady current may also vary upon dif- 
ferent time-dependent applied voltages, even if the volt- 
age amplitudes are same at i — J- -l-oo. For instance, there 
may be multiple steady states corresponding to the mul- 
tiple stationary solutions of the TDDFT-EOM, Eq. ([6]). 
The actual steady state reached at t — > +00 is thus com- 
pletely determined by the initial electron density, and the 
history of externally applied voltage. 

It has been shown that bound states in the device 
region would give rise to persistent oscillating current 
across the device-electrode interfaces4i^ In such cases, 
the Landauer-Biittiker formula becomes inadequate for 
describing the long-time limit of transient current. In 



fact, these bound states are physically decoupled with 
the rest of composite system, and form an isolated sub- 
system. When this isolated subsystem possesses nontriv- 
ial component (such as tail of wavefunction) inside the 
r-space of electrodes, the persistent oscillating current 
would arise, as an electron transfers across the device- 
electrode interfaces, but remains within the subsystem. 
It is thus inferred that the magnitude of oscillating cur- 
rent would diminish as the device region is enlarged by 
pushing the interfaces (where the currents are measured) 
deeper into the electrodes. The oscillating current would 
vanish completely with all the bound states fully accom- 
modated in the device region. 
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